浮点数的表示:详见二进制转换成其它格式数据,此处的浮点数精度为64位,即指数部分为11位,尾数部分为52位。标准化的IEEE浮点数表示为:
对于64位浮点数,精度为
IEEE舍入最近法则
如果在二进制数右边的第53位是0,则向下舍去(在第52位后面截断)。
如果第53位是1,则向上进位(在第52位上加1);但是如果在53位的右边的所有已知位都是0,此时当且仅当第52位是1时才在第52位上加1,否则在52位后面截断。
定义:将IEEE双精度浮点数字记做x,利用舍入最近法则记做fl(x)。
举例:对于数字9.4,二进制表示为
将二进制转成浮点数为:
对于
所以舍入误差为
相对舍入误差
在IEEE机器算术模型中,fl(x)的相对舍入误差不会比机器精度的一半大:
举例:fl(9.4)对应的存储在机器中的表达?
对于指数:3 + 1023 = 1026 =
对于尾数:直接转成16进制可得
整体16进制表示为4022CCCCCCCCCCCD。
特殊数的十六进制表示
| 机器数 | 例子 | 十六进制数 |
|---|---|---|
| +Inf | 1/0 | 7FF0000000000000 |
| -Inf | -1/0 | FFF0000000000000 |
| NaN | 0/0 | FFFxxxxxxxxxxxxx |
| +0 | 0000000000000000 | |
| -0 | 8000000000000000 |
其中x表示不全为0的位,+0和-0在计算中相同,只有符号位不同。
当指数全为0时,如0000000000000001表示的数为
对于1和
举例:计算双精度浮点的求和(

由于求和后的53位是1,又52位也是1,所以需要在第52位上加上1,即最后结果为:

在减去1后,结果为
举例:找到二次等式
使用求根公式:
这样便可解得正号的根
二分法
给定初始区间[a,b] 使得f(a)f(b)<0 while (b-a)/2 > TOL c = (a+b) / 2 if f(c) = 0, stop,end if f(a)f(c) < 0 b=c else a=c end end 最终的区间[a,b]中包含一个根 近似根为(a+b)/2
二分法的求解误差
其中
定义:如果误差小于
举例:使用二分法求解函数在区间[0,1]上的解,要求解精确到小数点后6位,则至少要迭代多少次。
至少要20次迭代。
定义:当g(r)=r,实数r是函数g的不动点。
不动点迭代
=初始估计
,其中i=0,1,2,...
举例:使用不动点迭代求解
将方程重写为
还可以将两边加上
使用
使用cobweb图描述不动点迭代的几何图示

可以看到对于g(x)=1-x^3而言,在迭代时会出现发散的情况,因此最后无法收敛,而是否能够收敛与函数在不动点附近的斜率有关。
定义:令
该方法被称为满足线性收敛,收敛速度为S。
定理:假设函数g是连续可微函数,g(r)=r,S=
证明:令
代入
定义
如果S=|g'(r)|小于1,则通过替换g',在r附近有一个小的区间满足|g'(x)|<(S+1)/2,这个值仍然小于1,如果
所以,误差以(S+1)/2的速度下降,这意味着
定义:如果迭代方法对于一个足够接近r的初值能收敛到r,该迭代方法被称为局部收敛到r。
举例:解释为什么不动点迭代g(x)=cos(x)收敛
g'(r)=-sin(r)=-sin(0.74)≈-0.67,绝对值比1小,对应解r≈0.74会把附近的初始猜测值吸引过来。
举例:使用不动点迭代找出方程 cos x = sin x 的根。
在方程两边加上x可得
定义g(x)=x + cos x - sin x,可用不动点迭代法求得x=0.7853981633957986。
举例:使用二分法计算函数
在使用二分法时,在16步计算后便出现错误,f(0.6666641)=0,无法精确到小数点后6位。即使使用MATLAB提供的fzero.m函数计算仍然无法找到5位以上的正确数位。
定义:假设f是一个函数,r是一个根,意味着满足f(r)=0。假设
后向误差在左边或者输入一侧(问题数据),这是我们需要对于问题(函数f)改变的量使得方程平衡并输出近似解
在上面的例子中,后向误差接近
方程求解方法的终止条件可以基于前向或者后向误差。还有其他相关的终止条件,诸如计算时间的上限问题的上下文对于我们的选择具有指导作用。
由于函数在多重根位置上f'为0,所以在重根附近形状很平。正因为如此,分离重根可能会遇到困难,但是多重根问题仅仅是冰山一角,没有多重根问题时也可能也会出现问题,比如威尔金森多项式。
如果在输人中是一个小误差,f在这种情况下对问题进行求解,造成输出中的大问题,这样的问题被称为敏感性问题。
假设对输入做了一个小变化εg(x),其中ε很小,令Δr对应根中的变化,则
进行泰勒展开得到:
忽略高阶项,得到
或者,由于ε很小,所以εg'(r)很小。
根的敏感公式
假设r是函数f(x)的根,并且r+Δr是f(x)+εg(r)的根,则当ε<<f'(r)时,
举例:估计
令
误差放大因子 = 相对前向误差 / 相对后向误差
条件数也是误差放大度量的一种方式,数值分析是对算法的研究,算法把定义问题的数据作为输入,对应的结果作为输出。条件数指的是理论问题本身所带来的误差放大部分,和用于求解问题的特定算法无关。注意到条件数仅仅度量由于问题本身带来的误差放大,这点很重要。和条件一起,还有一个平行概念,即稳定。稳定指的是由于算法小的输入误差造成的放大,而不是问题本身。如果一个算法在小的后向误差存在的时候,总能给出一个近似解,则称该算法是稳定的。如果问题的条件好,算法稳定,我们可以期望同时具有小的后向误差和前向误差。
前面的误差放大例子表明根求解过程对于特定的输入变化的敏感性。问题可能或多或少地敏感,依赖于如何设计输入的变化。问题的条件数定义为所有输入变化,或者至少规定类型的变化所造成的最大误差放大。条件数高的问题称为病态问题,条件数在1附近的问题称为良态问题。
牛顿方法
=初始估计
定义:令
定理:令f是二阶连续可微函数,f(r)=0。如果f'(r)≠0,则牛顿方法局部二次收敛到r,第i步的误差
证明参见原书,注意此处
注意,该定理并不能保证牛顿方法每次都能二次收敛。如求f(x)=
误差公式是
定理:假设在区间[a,b]上,(m+1)阶连续可微函数f在r点有一个m阶多重根,则牛顿方法局部收敛到r,第i步误差
定理:如果在[a,b]区间上f是(m+1)阶连续函数,包含m>1的多重根,则改进的牛顿方法为:
收敛到r,并具有二次收敛速度。
牛顿方法如同FPI(不动点迭代),可能不会收敛到根.下面的例子是一种可能不收敛的情况。
举例:对函数f(x)=
除了重根,牛顿方法比二分法和FPI方法的收敛速度更快。它达到了这种更快的速度是因为使用了更多的信息,尤其是通过函数导数得到的函数切线方向的信息。在某些情况下,可能难以计算导数。在这种情况下,割线方法是牛顿方法的一个非常好的替代。它使用近似值割线替代了切线,并且收敛速度差不多快。
在
割线方法
, =初始估计
割线方法以超线性的速度收敛到一个单根,意味着它在线性和二次收敛方法之间。
试位方法和二分法相似,将中点从a和b的平均值改为割线方法定义的点,注意f(a)f(b)<0
c对应的便是a和b所定义的直线与x轴的交点。
试位方法开始表现得比二分法和割线方法都要好,具有二者最好的性质。但是,二分法在每一步中可以确保消除1/2的不确定性,试位方法却没有能力做出这样的保证,而且在一些例子中可能收敛很慢。
Muller方法是割线方法在不同方向的推广。该方法不是计算经过先前两个点的直线和x轴的交点,而是使用三个前面生成的点
逆二次插值(IQI)是割线方法到抛物线的一种相近的泛化方法,但是使用形如x=p(y)的抛物线,而不是 Muller方法所使用的y=p(x),这个抛物线和x轴只有一个交点。经过三点(a,A),(b,B),(c,C)的二阶多项式x=P(y)为:
其中P(A)=a,P(B)=b,P(C)=c,用y=0带入可得得到抛物线和x轴的交点,可以得到
其中q=f(a)/f(b),r=f(c)/f(b),s=f(c)/f(a)。
逆二次插值
, , , , , 下一步的
为
可以用新的
替代最旧的估计 ,另一种方式使用最新估计替换最近的三个估计中的后向误差最大的那个。
用于连续函数f,区间的边界是a和b,同时f(a)f(b)<0。Brent方法记录当前点
(1)后向误差得到改进;
(2)包含根的区间至少减小一半;
否则,尝试使用割线方法以实现相同的目的。如果割线方法也失败了,则使用二分法,保证至少减少一半的不确定性。
高斯消去法不过多介绍
可以将高斯消去法的消去步骤表示为LU分解,其中L是单位矩阵进行高斯消去操作后的下三角矩阵,U是原矩阵经过高斯消去得到的矩阵。
一旦知道L和U,问题Ax=b可以写作LUx=b,定义新的“辅助”向量c=Ux,则回代是个两步的过程
(1)对于方程Lc=b,求解c
(2)对于方程Ux=c,求解x
但并不是所有的矩阵都能进行LU分解,如
无法被LU分解。
在高斯消去法中有两个潜在的误差来源,第一个是病态问题的概念和解对于输入数据的敏感性相关,第二个是淹没,在大部分问题中可以通过一个简单的修正,称为部分主元,进行避免。
定义:向量
定义:令
举例:找出近似解[-1,3.0001]的后向误差和前向误差,方程组如下:
精确解为[1,1],后向误差计算:
后向误差为0.0001,前向误差计算:
前向误差为2.0001。可以看到前向误差和后向误差差别很大。下面这张图可以帮助理解这一情况

可以看到(1,1)是它们的交点,不过(-1,3.0001)也差一点在两条直线上(图中的差异是人为放大了1000倍,实际上离得非常近)。
把余项表示为
相对前向误差定义为
误差放大因子则是二者的比值,即 相对前向误差/相对后向误差。
定义:方阵A的条件数cond(A)为求解Ax=b时,对于所有右侧向量b,可能出现的最大误差放大因子。
方阵A的矩阵范数定义为每行元素绝对值求和,其中的最大值作为矩阵A的范数。
定理:n×n矩阵A的条件数是:
举例:计算矩阵A的条件数,其中
计算
在浮点算术中,相对后向误差不可能小于
以上例为例,
对于希尔伯特矩阵H(
对于向量和矩阵范数,有
对于向量x的1-范数为
举例:考虑方程组
(1)精确解约为[2,1]
(2)IEEE双精度,
由于存在舍入,
(3)IEEE双精度,使用行交换
同样存在舍入导致
第2种方法在消去过程中的乘子
而另一方面,第3种方法在消去过程中没有产生覆盖,因为乘子是
高斯消去的过程中保证乘子尽可能小,同时避免淹没。幸运的是,通过对朴素的高斯消去法的简单修正,可以使得高斯消去法中的乘子的绝对值不大于1。这种新的原则,包括在表格中明智的行交换,称为部分主元方法,在下一节中将进行详细的描述。
部分主元包含在每一步消去步骤之前的比较,找到第一列中最大的一个元素,其对应行和主元行进行交换,在当前情况下,主元行是第一行。在消去过程中,对每一列都实施相同的策略。在消去第k列时,找到第p行,k≤p≤n,定位最大的
举例:使用高斯消去的部分主元方法求解方程组
由于第一列中
在消去第二列时,由于
解得x=[1,1,-1]
定义:置换矩阵是一个n×n的矩阵,其在每一行、每一列仅有一个1,其他全部为0。
定理:令P是通过对单位矩阵实施一组特定的行交换后得到的一个n×n的置换矩阵。则对于任意的n×n矩阵 A,PA对应于对矩阵A实施同样的行交换得到的结果。
我们把关于高斯消去所知道的一切组合起来,进行PA=LU分解。这是部分主元的高斯消去的矩阵形式,PA=LU分解是求解线性方程组的主要方法。正如名字中所暗示的,PA=LU是对于矩阵A包含行交换的LU分解。在部分主元中,开始的时候我们并不知道需要进行置换的列,所以我们必须谨慎地把行置换的信息加 入分解中,我们需要记录高斯消去法中所有的主元。
使用PA=LU求解方程组Ax=b时
其中P即为单位阵进行部分主元法交换行后产生的置换矩阵。
雅可比(Jacobi)方法是方程组系统中的一种形式的不动点迭代。在FPI中,第一步是重写方程,进而求解未知量。雅可比方法以如下标准方式进行该重写步骤:求解第i个方程得到第i个未知变量,然后使用不动点迭代,从初始估计开始,进行迭代。
举例:使用雅可比方法求解下面的线性方程组
求解未知变量得到
迭代过程为:
定义:当对角线上的元素的绝对值比该行上其它元素的绝对值之和还要大的时候,认为n×n的矩阵
定理:如果n×n矩阵A是严格对角占优矩阵,则(1) A是非奇异矩阵,(2)对于所有向量b和初始估计,对Ax=b应用雅可比方法都会收敛到(唯一)解。
但是这个只是充分条件,不代表不满足严格对角占优矩阵,就不会收敛到唯一解,上例中就不是严格对角占优矩阵,可以看到第一行的对角线元素为3,只是等于第一行其它元素的绝对值之和。
雅可比方法是不动点迭代的一种形式,令D表示A的主对角线矩阵,L表示矩阵A的下三角矩阵(主对角线以下的元素),U表示上三角矩阵(主对角线以上的元素),则A=L+D+U,求解的方程变为Lx+Dx+Ux=b。 注意,这里对于L和U的使用和LU分解中不同,当前L和U的主对角线元素都是零。
Ax=b
(D+L+U)x=b
Dx=b-(L+U)x
x=
令g(x)=
雅可比方法
=初始向量
= (b-(L+U) ), k = 0,1,2,...
高斯-塞德尔方法和雅可比方法的唯一差异是在前者中,最近更新的未知变量的值在每一步中都使用,即 使更新发生在当前步骤。在使用相同步骤的情况下,该近似更加精确。如果收敛,高斯-塞德尔方法常常比雅可比方法收敛更快。
高斯-塞德尔方法
=初始向量
= (b-U -L ) k = 0,1,2,...
举例:对如下系统应用高斯-塞德尔方法
迭代公式为
连续过松弛(SOR)方法使用高斯-塞德尔的求解方向,并使用过松弛以加快收敛速度。令w是一个实数,并将新的估计中的每个元素
使用SOR后的迭代公式
可以将系统看成不动点迭代,如下:
连续过松弛(SOR)
=初始向量
收敛性的证明略
基于高斯消去的直接方法,通过有限步数的计算就可以得到解。为什么我们要研究并使用迭代方法,特别是迭代方法只能得到近似解,并需要多步计算以得到收敛。有两个主要的原因让我们使用诸如高斯-塞德尔这样的迭代方法。两个原因都是基于如下的事实:迭代方法中的一步计算,仅仅需要LU分解的浮点计算所需要时间的一小部分。对于n×n的矩阵做一次高斯消去需要
如果已知解的较好的近似,在这种特定的情况下支持使用迭代技术。例如,假设知道Ax=b的解,随后A与b同时或者单个仅仅发生小的变化。我们可以想象一个动态系统,当A和b改变时,对A和b进行待续的度量,因而需要待续更新精确的解向量x。如果把前面问题的解作为新的相似问题的初始估计,使用雅可 比或者高斯-塞德尔可以得到更快的收敛。这种技术通常称为修饰,这是由于该问题从一个近似解开始,该近似解可能对应前面一个相关问题的解,然后仅仅修饰近似解使其更加精确。
第二个使用迭代方法的主要原因是用于求解稀疏的方程组。当矩阵中的很多元素都是0系数矩阵被称为稀 疏矩阵。通常,对于稀疏矩阵中的
x% 程序 2.1 稀疏矩阵生成% 输入:n = 系统大小% 输出:稀疏矩阵 a, r.h.s. bfunction [a,b] = sparsesetup (n)e = ones(n,l); n2=n/2;a= spdiags( [-e 3*e -e], -1 : 1,n,n); % a的元素c= spdiags( e/2, 0,n,n) ; c= fliplr(c); a= a+c;a(n2+1,n2) = -1; a (n2,n2+1) = -1; % 设置2个元素b=zeros(n,1); % r.h.s. b的元素b(1)=2.5; b(n)=2.5; b(2:n-1)=1.5; b(n2:n2+1)=1;
% 程序 2.2 雅可比方法% 输入:完全矩阵或者稀疏矩阵 a, r.h.s. b,% 雅可比迭代次数k% 输出:解xfunction x = jacobi(a, b, k)n=length(b); % 确定nd=diag(a); % 提取a的对角线元素r=a-diag(d ); % r为余项x=zeros(n,1); % 初始化向量xfor j =1:k % 雅可比迭代的循环 x = (b-r*x) ./ d;end
定义:如果
性质:如果n×n矩阵A是对称矩阵,则A是正定矩阵,当且仅当所有特征值是正数。
性质:如果A是n×n对称正定矩阵,X是一个满秩n×m矩阵,n≥m,则
定义:方阵A的主子矩阵是一个方的子矩阵,其对角线元素是矩阵A的对角线元素。
性质:任何对称正定矩阵的主子矩阵也是对称正定矩阵。
考虑一个对称正定矩阵
将A写成
与上式作比较,
定理:如果A是n×n对称正定矩阵,则存在上三角n×n矩阵R满足
Cholesky分解
for k = 1,2,...,n
if
<0,stop,end
=
end
举例:计算如下矩阵的Cholesky分解
k=1:
k=2:现在在这个2×2的子矩阵中重复这个过程,找到
共扼梯度的思路依赖于内积思想的推广,因为(v,w)=(w,v),以及对于标量α和β,有(αv+βw,u)=α(v,u)+β(w,u),所以欧几里得内积
定义:令A是对称正定的n×n矩阵,对于两个n维向量v和w,定义A内积
当
共轭梯度方法
=初始估计
for k = 0,1,2,...,n-1
if
=0,stop,end
end
对于
对于
有
通过使用预条件技术,可以使得诸如共扼梯度方法的迭代方法收敛速度加快。迭代方法的收敛率通常直接或者间接依赖于系数矩阵A的条件数。预条件方法的思想是降低问题中的条件数。
n×n的线性方程组Ax=b的预条件形式是
其中M是可逆的n×n矩阵,称为预条件子。我们所要做的就是在方程两侧左乘上该矩阵,预条件矩阵试图对矩阵A逆转,从而可以有效地降低问题的条件数。概念上来讲,它试图同时做两件事:矩阵M应该(1)和矩阵A足够接近,(2)容易求逆。这两个目的常常彼此对立。
一种特别简单的方式是雅可比预条件子M=D,其中D是A的对角线矩阵。D的逆矩阵是D的元素的倒数.例如在一个严格对角占优矩阵中,雅可比预条件子和A相似,同时非常容易求逆。注意到,对称正定矩阵的每个对角线元素都严格为正,因而计算倒数不是问题。
预条件共轭梯度法
=初始估计
for k = 0,1,2,...,n-1
if
=0,stop,end
end
在对称连续过松弛(SSOR)中,预条件子定义如下:
其中A=L+D+U被分为下三角部分、对角线以及上三角部分。w是一个在0和2之间的常数。在特例中w=1,这被称为高斯-塞德尔预条件子。
注意到SSOR预条件子定义为下三角矩阵和上三角矩阵的乘积
多变量情况下函数导数f'对应的是雅可比矩阵,定义为
多变量牛顿方法
=初始向量
,k=0,1,2,...
由于出现计算矩阵的逆,所以令
如果可以计算雅可比矩阵,那么牛顿方法是一个好选择,那么最好的替代方法就是Broyden方法。
Broyden方法
=初始向量
=初始矩阵 for i = 0,1,2,...
end
其中
Broyden方法Ⅱ
=初始向量
=初始矩阵 for i = 0,1,2,...
end
其中
首先,需要初始向量
定义:如果对于每个1≤i≤n,
假设给出三点(
一般来说,如果有n个点
定理(多项式插值的主定理):令
定义:用f[
牛顿差商公式:
牛顿差商
给定x=[
,..., ],y=[ ,..., ] for j = 1,...,n
f[
]= end
for i = 2,...,n
for j=1,...,n+1-i
f[
... ]=(f[ ... ]-f[ ... ])/( - ) end
end
xxxxxxxxxx% 程序 3.1 牛顿差商插值方法% 计算插值多项式的系数% 输入:x 和y是包含 n 个数据点的 x 和 y坐 标的向量% 输出:嵌套形式的插值多项式系数 c% 使用 nes 七 . m 计算插值 多项式function c=newtdd(x, y,n)for j =1:n v(j,1) =y (j) ; % 填入牛顿三角形的y列endfor i =2:n % 对于第l.列 for j =1 : n+1- i % 从顶端到底端填入列元素 v (j, i) = (v (j + 1, i -1) - v(j, i -1)) / (x (j +i - 1) -x (j)); endendfor i = 1 : n c(i) =v(1, i) ; % 从三角形顶端读end % 输出系数
不过多赘述。
切比雪夫插值是一种特定最优的点间距选取方式。
切比雪夫插值的动机是在插值区间上,提高对如下插值误差的最大值的控制
从现在开始让我们把区间固定在[-1,1]。多项式插值误差的分子
本身是一个关于x的n阶多项式,并在区间[-1,1]上具有极值。是否可能在区间[-1,1]找到特定的
定理:选择实数
尽可能小,则
对应的最小值是
可以得到极小值,其中
定义n阶切比雪夫多项式
事实:
事实:
事实:
事实:
事实:
关于切比雪夫插值的讨论局限在区间[-1,1],可以将方法推广到一般的区间[a,b],移动基点使得它们在区间[a,b]上的相对位置和在区间[-1,1]上一致。这可以通过如下两步实现:
(1)使用因子(b-a)/2拉伸点(这是两个区间长度的比值),
(2)将点平移(b+a)/2,使得中心从0移动到区间[a,b]的中心。
换句话讲,从原始点
使用区间[a,b]上新的切比雪夫基点
切比雪夫插值节点
在区间[a,b]
i=1,...,n,不等式
在区间[a,b]上成立。
举例:在区间[0,π/2]找出4个切比雪夫基点进行插值,找出切比雪夫插值误差的上界,在区间中f(x)=sinx。
插值误差的最坏情况是
这里获得的只是插值时使用的基点,即牛顿差商中
样条是另一种数据插值的方式,在多项式插值中,多项式给出的单一公式满足所有数据点,而样条 则使用多个公式,其中每个都是低阶多项式,来通过所有数据点。线性样条可以成功地对任意的n 个点集进行插值,但是线性样条函数缺乏平滑,三次样条则可以解决线性样条的这个缺点。三次样条在两个数据点之间使用3阶(cubic)多项式替换线性样条两点之间的线性函数。由于通过n个点的三次样条无穷多,所以需要添加额外条件。
假设给定n个点
并具有如下性质:
(1)
(2)
(3)
性质1保证样条S(x)插值数据点;性质2使得相邻的样条段在它们相遇的地方斜率相同;性质3则保证在两条样条段相邻的地方曲率相同,该曲率由二阶导数表示。
只有满足上面的这些条件才被称为三次样条
性质4a(自然样条):
满足性质4a这两个附加条件的三次样条被称为自然三次样条,自然三次样条是唯一的。

自然三次样条
给定
,其中 , for i = 1,...,n-1
end
求解3.24得到
for i = 1,...,n-1
end
除了性质4a给出的自然三次样条,还有曲率调整三次样条、钳制三次样条、抛物线端点的三次样条、非纽结三次样条。
贝塞尔样条是一个允许用户控制节点处斜率的样条。作为额外自由控制的代价,在节点处的一阶和二阶导数的平滑性不再能保证,而这种平滑性是前面三次样条本身就具有的性质。贝塞尔样条适合不时出现角点(一阶导数不连续)和曲率突变(二阶导数不连续)的情况。
平面贝塞尔样条的每一段由4个点

贝塞尔曲线
给定端点
, 控制点
, 设
定义在0≤t≤1的贝塞尔曲线如下:
计算机中的字母或者数字便是使用二维贝塞尔曲线画出来的,如Times Roman字体中大写T字符由下面的16条贝塞尔曲线构成,每条曲线包含
237 620 237 620 237 120 237 120; 237 120 237 35 226 24 143 19; 143 19 143 19 143 0 143 0; 143 0 143 0 435 0 435 0; 435 0 435 0 435 19 435 19; 435 19 353 23 339 36 339 109; 339 109 339 108 339 620 339 620; 339 620 339 620 393 620 393 620; 393 620 507 620 529 602 552 492; 552 492 552 492 576 492 576 492; 576 492 576 492 570 662 570 662; 570 662 570 662 6 662 6 662; 6 662 6 662 0 492 0 492; 0 492 0 492 24 492 24 492; 24 492 48 602 71 620 183 620; 183 620 183 620 237 620 237 620;
画出来的曲线为:

python代码: BezierCurve[贝塞尔曲线].py
在方程组一章中,给出了当方程解存在时,如何找到Ax=b的解。我们将学习当解不存在的时候该怎么办。当方程不一致时,有可能方程的个数超过未知变量的个数,答案是找到第二可能好的解:即最小二乘近似.
在插值一章中给出了如何找出多项式,并精确拟合数据点。但是如果有大量的数据点,或者采集的数据点具有一定误差,使用高阶多项式精确拟合一般不是个好方法。在这种情况下,使用简单模型 近似数据是一种更合理的方式。
求解不一致的方程组以及近似拟合数据这两个问题一同驱动着最小二乘的发展.
如果一个方程组无解,它被称为不一致。方程无解意味着什么?可能系数不十分精确,在很多情况下,方程的个数比未知变量的个数要多,使得解不可能满足所有的方程。但是还可以通过找出与解最相似的向量x。
如果我们选择“相似度”意味着欧氏距离相近,有一个直接的方法找到最接近的x,这个特殊的x称为 最小二乘解。
对于Ax=b这个方程而言,如果b不在Ax可以确定的平面(A中的每一列看成基向量,对应的x的一行(就是一个值)便可看成该方向的值)上,那么就是无解的,不过可以在这个平面上找到与b最接近的向量来近似。
最小二乘基于正交,从一点到一个平面的最短距离,由一个到平面的正交线段表示,法线方程可以确定该线段,这表示最小二乘的误差。
假设这个近似解为
此时这个方程便是可解的,不过解出的结果与实际结果可能存在一定差距,有三种方式可以表示余项的大小。
(1)欧式长度:
(2)2范数,平方误差:
(3)平均平方根误差:
三种表达紧密相关。
最小二乘数据拟合
给定一组数据点
,..., step1 选择模型,确定用于拟合数据的参数模型,例如y=a+bt
step2 将数据点代入模型,每个数据点生成一个方程,其中的未知变量是模型参数,例如线性模型中的
和 为参数,这得到系统Ax=b,其中未知变量x表示未知参数。 step3 参数的最小二乘解是法线方程
的解。
最小二乘是数据压缩的经典例子,输入包含一组数据点,输出是一个尽可能好的数据拟合模型,其中具有较少的参数。通常,使用最小二乘的原因是使用合理的底层模型替换噪声数据,该模型通常用于信号预测以及分类。
当
周期数据使用周期模型。例如外层大气温度遵循大时间尺度的循环,包含由地球自转和绕太阳公转控制的每天和每年的循环。
如拟合一天为周期的温度变化,可以使用模型
以描述人口变化的指数模型为例:
不能直接进行最小二乘拟合,这是由于
不过有两种方法可以用来解决这个问题,第一种方法是直接求解非线性最小二乘方程(在后面介绍),第二种是通过取对数来使问题线性化。
同样对于幂法则模型
格拉姆-施密特方法是对一组向量正交化。给定一组输入的m维向量,目的是找出正交坐标系统,获取由这些向量张成的空间。更精确地讲,给定n个线性无关的输入向量,该方法计算n个彼此垂直的单位向量(单位长度由2范数定义),张成和输入向量相同的子空间。
格拉姆-施密特正交的思想是对向量减去其投影在其它已经正交好的向量上的分量,得到正交于这些正交好的向量,在进行单位化即可。令
或者写成矩阵形式
或A=QR,其中A是包含列向量
经典格拉姆-施密特正交
令
(j=1,...,n)为线性无关向量。 for j = 1,2,...,n
for i = 1,2,...,j-1
end
end
当成功分解后,通常会填满正交单位向量的矩阵,从而得到
注意A是m×n矩阵,Q是m×m方阵,上三角矩阵R是m×n矩阵,m>n,这称为A的完全QR分解,Q在数值分析中具有特殊的地位,并具有一个特殊的定义。
定义:当
引理:如果Q是m×m正交矩阵,x是m维向量,则
LU分解是对高斯消去中信息进行有效编码的方式。以相同方式,QR分解记录了矩阵正交化的信息 ,即构造一个正交集,张出由A的列向量构成的空间。使用正交矩阵计算更好的原因是
(1)根据定义它们很容易求逆,
(2)由引理知,它们不会放大误差。
QR分解可用于求解n个方程n个未知量的方程组,但是计算复杂度比LU分解大。
QR分解可以用于最小二乘,令A是m×n矩阵,其中m≥n。为了最小化
欧几里得范数中的向量是
其中
通过QR分解实现最小二乘
给定m×n不一致系统
找出完全QR分解A=QR,令
的上n×n子矩阵
的上面的n个元素 求解
得到最小二乘解
相比于法线方程求解最小二乘,这个更精确。,可以更好处理病态问题。
改进的格拉姆-施密特正交
令
(j=1,...,n)为线性无关向量。 for j = 1,2,...,n
for i = 1,2,...,j-1
end
end
与经典格拉姆-施密特唯一不同的是,
尽管改进的格拉姆-施密特正交是计算矩阵的QR分解的有效方式,但是这不是最好方式。另一种方法是使用豪斯霍尔德反射,这种方式需要更少的计算,同时在舍入误差放大的意义上来讲这种方法也更稳定。
豪斯霍尔德反射子是正交矩阵,通过m-1维平面反射m维向量。这意味着当每个向量乘上矩阵后,长度保持不变,使得豪斯霍尔德反射被称为移动向量的完美形式。给定一个向量x,我们要重新找出一个相同长度的向量w,计算豪斯霍尔德反射得出矩阵H满足Hx=w。
如图所示,原始方法非常清晰。画出m-1维平面二分x和w,并和连接它们的向量垂直,然后通过该平面反射所有向量。

引理:假设x和w是具有相同欧几里得长度的向量,
定义向量v=w-x,考虑投影矩阵
投影矩阵是一个满足
矩阵H被称为豪斯霍尔德反射子,注意到H是对称并正交的矩阵。
定理(豪斯霍尔德反射子):令x和w是向量,
使用豪斯霍尔德反射子推导出QR分解的一个新方式,使用反射子的目的是将列向量x移动到坐标轴,并以此将0放在矩阵中。
我们从矩阵A开始,希望写做形如A=QR的方程。令
我们已经在A中引入一些0,我们希望继续这种方式直到A变为上三角;然后我们将得到QR分解中的R。对于一个3列的矩阵,还需要做两次这样的操作,即
举例:使用豪斯霍尔德反射子找到矩阵A的QR分解
豪斯霍尔德反射子
下面使用
因此,有
与格拉姆-施密特正交相比,QR分解的计算代价更低,可以更好地得到单位正交向量以及具有更低的内存需求。
GMRES属于Krylov方法,这些方法依赖精确的Krylov空间的计算,该空间是向量{r,Ar,...,
GMRES的思想是在特殊矢量空间,即Krylov空间寻找初始估计的
广义最小余项方法
=初始估计
for k = 1,2,...,m
for j = 1,2,...,k
end
(如果 ,跳过下一行,并在底端终止)
最小化
得到
end
迭代的
在预条件GMRES背后的概念和共枙梯度法非常相似,从非对称的线性方程组Ax=b开始,我们试图求解
预条件GMRES
=初始估计
for k = 1,2,...,m
for j = 1,2,...,k
end
最小化
得到
end
首先会议多变量的牛顿方法,并将其应用在列向量函数
其中
高斯牛顿方法
为了最小化
令
=初始向量, for k = 0,1,2,...
end
定义矩阵
举例:拟合模型
定义余项
并相对
最小二乘问题中的非线性带来额外的挑战。法线方程以及QR方法只要系数矩阵A满秩都可以找到唯一解。而对于非线性问题的高斯-牛顿迭代可能收敛到多个极小平方误差中的一个。尽可能使用初始向量的合理近似,有助于收敛到绝对极小。
非线性最小二乘法遇到病态系数矩阵时,会面临巨大的挑战,Levenberg-Marquardt方法使用“正则化项”部分修复这个问题。
Levenberg-Marquardt方法
最小化
令
=初始向量, =常数 for k = 0,1,2,...
end
提高正则化参数λ可以增强矩阵
计算微积分解决的主要问题是计算函数的导数和积分。
二点前向差分公式
其中c在x和x+h之间,
定理(推广中值定理):令f是区间[a,b]上的连续函数,令
三点中心差分公式
其中x-h<c<x+h。
二阶导数的三点中心差分公式
其中x-h<c<x+h。
当h->0时,误差会先减小再增大,后增大的原因是计算时有效数字的丢失。对于三点中心差分公式而言,f'(x)的及其近似误差的绝对值上界是
假设我们有n阶公式F(h)近似一个给定量Q,这个阶数意味着
其中K大约是我们感兴趣的h区间上的一个常数,相关的例子是
n阶公式的外推
举例:应用外推公式
梯形法则
其中
辛普森法则
其中
定义:数值积分方法的精度是最大的整数k,使用该积分方法可以得到所有k阶或者更低阶多项式积分的精确值。
梯形法则的精度是1,辛普森法则的精度是3。
由于积分在区间的所有子区间上具有可加性,我们可以通过除法把整个区间变为很多小区间再计算积分,在每个小区间上使用法则,然后再求和。这种策略被称为复合数值积分。
复合梯形法则
其中h=(b-a)/m,c在a和b之间。
复合辛普森公式
其中c在a和b之间。
中点法则
其中
中点法则也可用于减少所需函数求值的次数与梯形法则相比,闭牛顿-科特斯方法具有相同的阶数,只需要一次函数求值,而不是两次。而且误差项是梯形法则误差项大小的一半。
复合中点法则
其中h=(b-a)/m,c在a和b之间,
龙贝格积分是对复合梯形法则应用外推的结果,给定近似M的法则N(h),该法则依赖步长h,如果知道法则的阶,则可以对法则进行外推。
龙贝格积分
for j = 2,3,...n
for k = 2,...,j
end
end
最后的结果便是
到目前我们学到的积分近似方法都使用相等的步长。一般来说,小步长可以提高精度变化剧烈的函数将需要更多步,因而带来更多的计算时间,这是由于需要更小的步长跟踪函数所有的变化。函数可能在定义域某些部分变化剧烈,而在另外一些部分变化缓慢,满足前一部分的容差的步长在后一部分就显得过于精细了。
通过使用积分误差公式,可以在运算中推出一个标准,其步长对于特定的子空间适合。这种方法称为自适应积分。
自适应积分
对给定容差TOL近似
c = (a+b)/2
=(b-a)(f(a)+f(b))/2 if
接受
作为区间[a,b]上的近似 else
对于区间[a,c]和[c,b]递归重复上面的步骤
end
xxxxxxxxxx% 程序 5 . 2 自适应积分% 计算定积分的近似% 输入: MATLAB函数f,区间 [ ao, bo]% 容差 tol0% 输出:近似定积分function int=adapquad(f, a0, b0, tol0)int = 0; n=1; a(1) = a0; b(1)=b0; tol(1)=tol0; app(1)=trap(f,a,b); % trap求梯形积分while n >0 % n是当前列表结束的位置。 c = (a(n) + b(n)) / 2; oldapp = app(n); app(n) = trap(f, a(n), c); app(n+1) = trap(f, c, b(n)); if abs(oldapp - (app(n) + app(n+1))) < 3 * tol(n) int = int + app(n) + app(n+1); % 成功 n = n-1; % 该区间操作完成 else % 分成两个子区间 b(n+1) = b(n); b(n) = c; % 设置新区间 a(n+1) = c; tol(n) = tol(n) / 2; tol(n+1) = tol(n); n = n + 1; % 到列表尾端,重复 endend其中tol0是自己设置的参数,如设置为0.005。
定义:一组在区间[a,b]上的非0函数{
定理:如果{
定理:如果{
高斯积分
其中
上面带横线的项表示不进行计算,
定理:高斯积分方法,在区间[—1,1]上使用n阶勒让德多项式,具有2n-1阶精度。
一般区间[a,b]的积分可以通过替代t=(2x-a-b)/(b-a)得到
通过跟随箭头计算“求解”微分方程。从初始条件(
欧拉方法
h是步长,如取0.2,y'=f(t,w)。
举例:用欧拉方法求解如下问题
欧拉方法则对应如下迭代,h设置为0.1
定义:当存在常数L(称为利普希茨常数)对矩形 S=[a,b]×[α,β]中的每对(t,
函数f(t,y)相对于变量y在S上利普西茨连续。
定理:假设f(t,y)定义在集合[a,b]×[α,β]并且α<
有唯一解y(t)。而且,如果f在[a,b]×(-∞,∞)是利普西茨连续,则它在[a,b]上存在唯一解。
定理:假设f(t,y)在集合S=[a,b]×[α,β]上关于y是连续的,如果Y(t)和Z(t)是微分方程
在S上的解,分别具有初值条件Y(z)和Z(a),则
一组容易求解的特定的常微分方程提供了ODE求解问题中具有指导意义的例子,这组方程是一阶方程,其右侧是关于y的线性函数。考虑初值问题
首先注意到如果g(t)在区间[a,b]上连续,唯一解存在,使用
积分因子是
显式梯形方法
泰勒方法基本思想是直接利用泰勒展开,假设解y(t)是(k+1)阶连续可微的函数。给定在解曲线上的当前点(t,y(t)),目标是对于某个步长h,使用微分方程的信息,用y(t)来表达y(t+h)。
K阶泰勒方法
举例:对于一阶线性方程确定二阶泰勒方法
由于
所以二阶泰勒方法的迭代公式为
单个的高阶微分方程可以转化为一个方程组,令
是n阶常微分方程,定义新的变量
注意到原始的常微分方程写成
两者放在一起
这些方程将n阶的微分方程转化为一阶微分方程组,该方程组可以使用欧拉或者梯形方法求解。
我们已经知道一阶的欧拉方法以及二阶的梯形方法。除了梯形方法外,还有其他龙格-库塔类型的二阶方法。
中点方法
4阶龙格-库塔方法(RK4)
其中
直到现在,步长h在ODE求解器实现中一直是一个常数。但是,并没有原因说h在求解过程中不能改变。我们希望改变步长的一个原因是,问题的解会在缓慢变化的周期和迅速变化的周期之间移动。将固定步长变得足够小使其精确跟踪快速变化,意味着解的其他部分求解都会是令人难以忍受的缓慢。
可变步长方法的关键思想是检测当前步生成的误差。用户设置的容差在当前步必须能够满足.然后设计方法(1)如果超出容差,必须拒绝误差并将步长减小,或者(2)如果满足容差,接受步长并选择对于下一步适合的步长h。关键是近似在每步中产生的误差。首先假设我们已经找到这样的方式并解释如何改变步长。
有RK2/3(龙格-库塔2阶/3阶嵌入对)、RK4/5,不过多赘述。
后向欧拉方法
举例:对初值问题使用后向欧拉方法
使用后向欧拉方法
令
可以使用牛顿方法进行估计
用
我们已经研究过的龙格-库塔方法家族包含单步方法,意味着新一步的
Adams-Bashforth两步方法
Adams-Bashforth三步方法(三阶)
Adams-Bashforth四步方法(四阶)
隐式梯形方法(二阶)
Adams-Moulton两步方法(三阶)
Milne-Simpson方法
Adams-Moulton三步方法(四阶)
Adams-Moulton四步方法(五阶)
令y(t)是一个至少四阶连续的函数,可以得到一阶导数的离散近似为
二阶导数的近似为
有限差分方法包含使用离散版本替换微分方程中的导数,并求解得到的更简单的代数方程获取正确值
举例:使用有限差分求解
使用二阶导数的中心差分形式,在
或者等价地
当n=3时,区间大小h=1/(n+1)=1/4,再插入边界条件
代入h得到三对角线矩阵方程
对于非线性边值问题,使用多变量牛顿方法进行迭代,迭代公式为
一般地,完成迭代的最好方式是求解方程
举例:求解非线性边值问题
离散化后方程为
其中
函数F(w)如下
和有限差分方法相似,排列和有限元方法背后的思想是将边值问题消减为一组可求解的代数方程。但是,并不是通过使用有限差分替换微分方程中的导数进行离散化,而是对于解给出函数形式,其对应的参数使用该方法拟合。
选择一组基函数
找出近似解的问题被简化为确定
考虑BVP问题
选择n个点,从边值点a开始,在b点结束,即
排列方法将的候选解代入微分方程中,并估计微分方程在n个点上的值,得到关于n个未知变量
选择基函数
所以可以写出n个未知变量
余下的n-2个方程来自微分方程在
在
在排列方法中选择三角函数作为基带来了傅里叶分析以及谱方法,这在边值方程和偏微分方程中都大量应用。这是一个“全局的方法,其中对于t的一个非常大的范围中,基函数非0,但是具有很好的正交性质。
选择样条函数作为基函数带来有限元方法.在这种方法中,每个基函数仅在t的一个小区间上非0。有限元方法大量用于高维情况下的BVP和PDE,特别是当存在不规则的边界条件的情况下,其中使用标准的基函数做参数化变得困难。Galerkin方法最小化微分方程在解上的平方误差。这带来关于
考虑一个BVP问题
解为y,余项r=y"-f,需要使得微分方程两侧的差异尽可能小。和最小二乘方法类似,这可以通过选择y使得余项与潜在解的向量空间正交得到。
对于区间[a,b],定义平方可积分的函数对应的向量空间
当
Galerkin方法包含两个主要思想,第一个是最小化r,其中强制其在
其中0≤i≤n+1,该形式被称为边值问题的弱形式。
Galerkin方法的第二个思想是部分使用积分消去二阶导数,注意到
可以得到
以如下函数形式求解
Galerkin方法的两个思想使其可以方便地使用极简单的函数作为有限元
同时定义
对于
同时有
第一个和最后一个
对于i=1,...,n,使用有限元方程
代入
定义偏微分方程(PDE)如下
使用数值微分和积分的离散公式近似在x和t方向的导数。例如,使用中心差分公式,关于x的二阶导数如下
误差为
其中误差
通过时间步进可以得到(定义
写成矩阵形式
此时离散的差分公式为(定义
矩阵形式为
Crank-Nicolson对于时间导数使用后向差分公式,并均匀组合前向差分近似和后向差分近似。
使用如下的后向差分公式替换
使用混合差分替换
可以重新组织热方程近似
或者
本章中有各种偏微分方程(抛物线方程,双曲线方程,椭圆方程,非线性偏微分方程)的解法。
定义:线性同余生成器(LCG)表示为以下形式(
其中a为乘子,b为偏移,m为模数。
最小标准随机数生成器
其中
, ,
按照现在计算机每秒指令数的增加,这个数目显得太少了。
randu生成器
其中
,
不满足随机数的独立性要求。
最新版的matlab已经不再使用LCG作为随机数生成器,从matlab5开始,使用的是延迟斐波那契生成器。
指数随机变量的取值满足概率分布函数
算法的核心思想是要找到指数随机变量的值,使得概率Prob(V≤x)是[0,1]区间上的均匀分布,即给定一个满足均匀分布的变量u之后,要使得
求解x得到
把输入的均匀随机数u转换为指数随机数。
对于一般的概率分布,也可以用这种方法处理。假设P(x)是需要生成的随机变量对应的概率分布函数。
标准正态分布也可以采用反函数的方式计算,不过还有一种更有效的方法,可以同时生成两个正态分布随机数。二维标准正态分布的概率密度函数为
其中
这种算法称为Box-Muller方法,此外还有一种更加高效的方法
该方法可以避免使用三角函数,但存在一定的拒绝率(21%左右),如果
Ⅰ型蒙特卡罗问题是用随机采样点计算方程的均值,再乘以积分区间的体积。计算函数均值可以看成是计算一个符合同样分布函数的概率的均值。我们用E(X)表示随机变量X的期望随机变量X的方差为
Ⅰ型或Ⅱ型蒙特卡罗模拟具有伪随机数字。
拟随机数的理念是在条件允许的情况下,放弃对随机数独立性的要求。放弃独立性意味着拟随机数不再是随机的,而且也不是像伪随机数那样看上去是随机的。通过这种方式,可以使得蒙特卡罗方法能很快地收敛到正确解。拟随机数的设计是要使得生成的序列是自避(Self-avoiding)的,而不一定要保证独立性。自避定义为,在生成随机数序列的过程中,新的数会填充到比较稀疏的区域,而不会聚集到一起。
Halton序列:设p为一个素数,例如p=2。按p进制表示写出1~n的自然数,假设第i个数字的p进制表示为
如p=3生成的前8个随机数
| i | i | ||||||
|---|---|---|---|---|---|---|---|
| 1 | 1 | 0.1 | 5 | 12 | 0.21 | ||
| 2 | 2 | 0.2 | 6 | 20 | 0.02 | ||
| 3 | 10 | 0.01 | 7 | 21 | 0.12 | ||
| 4 | 11 | 0.11 | 8 | 22 | 0.22 |
拟随机数的优势在于它的收敛速度更快。如果把估计误差写成计算次数n的函数的话,使用拟随机数会以更高阶的速度收敛。以下是采用拟随机数的误差收敛函数,d表示将要生成的随机数的维度:
1型蒙特卡洛问题(拟随机数)
2型蒙特卡洛问题(拟随机数)
随机游走模型,又称离散布朗运动。随机游走
其中t=0,1,2,...。
每个时刻t,
很多随机游走模型的应用都是关于逃逸时间,又称为首次到达时间。设a,b为正整数,一个初始位置为0的随机游走序列,首次达到[-b,a]区间边缘的时刻,就称为逃逸时间。理论证明在a处(而不是-b处)逃逸的概率为b/(a+b)。对于从区间[-b,a]逃逸所需的时间,其期望值为ab。
标准随机游走过程在时刻t的期望为0,方差为t。假设在每个单位时间内都可以进行两次位移,且每次位移的时间为1/2个单位时间,则在时刻t的期望仍然会是0,但是方差会变为
这是因为一共进行了2t次位移操作。如果我们把移动次数增加到原来的k倍,我们就必须把每次移动的距离缩小到原来的
因此,
当k趋向∞,极限
(1)对任意t,随机变量
(2)对任意
(3)布朗运动
定义:以实数t≥0为索引的随机变量集合
考虑如下的有噪声的微分方程(SDE)问题:
其中r和σ为常数,这个概率过程为
SDE是以微分的形式给出的,而在ODE中使用的是导数形式。这是因为像包括布朗运动在内的很多随机过程,它们是连续但不可微的。因此,SDE
等价于积分形式的表示
公式中第二个积分项称为Ito积分。
设
其中
其中
布朗运动
Ito公式
若y=f(t,x),则
有许多数值方法用于求解SDE,有Euier-Maruyania方法,Milstein方法,一阶随机龙格-库塔方法,不过多赘述。
DFT和快速傅里叶变换的内容不过多赘述。
设[c,d]为参数区间,n为正整数。定义
设y=
以上公式可以看成是对采样点
假设
把插值函数Q(t)=P(t)+I(t)分成实部和虚部两个部分,由于
下标n表示三角插值公式中的项数,有时可以把
当n为偶数时,有
满足
当n为奇数时,有
使用三角函数计算最小二乘近似问题
因为
定理:令
为实数正交阵,对y=Ax,以下函数
是经过采样点
令y=Ax,可以直接给出插值函数
DFT方法可以对区间[0,1]上的均匀采样点进行三角插值,假设n为偶数
公式中的项数为n,等于采样点的数目。采样点越多,需要用到的sin和cos函数项就越多。
当样本点数目 n 很大时,一般不会去精确地拟合模型函数.
实际上,在一般应用中,需要丢失一些信息(有损压缩)而使得模型比较简单;此外,由于采样点本身就是不准确的,因此强制要求插值函数经过所有点是不恰当的。对于这两种情况,可以使用上式进行最小二乘拟合,由于在这个模型中
令n表示样本点
当m<n时,就变成了一个压缩问题,我们希望用
最小二乘问题需要找到系数
的误差尽可能小,用矩阵形式表示,为
其中
n个点的插值多项式可以表示为
则解可以表示为
这个结论说明,对给定的n个数据点,用m<n个三角函数进行最小二乘拟合,只需计算n阶插值,再取前m个项即可。即x的插值系数Ax,在丢弃高频分量后,仍然是x的最佳拟合函数。n阶展开里的前m项,就是使用m个低阶频率可以达到的最佳拟合。这个性质反映了基函数的“正交性”。
推论:令[c,d]为一区间,m<n为正偶数,x=(
其中j=0,...,n-1,则有(注意
是数据
滤波的用法有两种。第一种是用另一个简单的函数尽可能地拟合原始音频.这是信号的压缩,我们可以只储存m个低频分量。滤波的另一个重要应用是去噪,在一个音乐文件里,音乐和语音中可能会混有高频噪声,需要消除这些高频分量来提高音频的质量。
设c为一个没有噪声的音频信号,与一个同样长度的向量r相加。得到的结果x=c+r是否是有噪声的 ?若r=c,我们认为r不是噪声,因为相加的结果只是增大c的音量,但信号是同样干净的。根据定
义,噪声是与信号无关的,即若r为噪声,则内积
在一个典型应用中,我们需要从含噪声的信号x中恢复出c。信号c可能是一个重要的系统变量,但是在噪声环境下进行测量。或者如下面例子,c可能是音频采样,我们想从中去除噪声。在20世纪中期Norbert Wiener建议,以最小平方误差为目标,寻找一种最优滤波器来去掉x中的噪声。他建议寻找一个对角阵Φ,使得公式
的模尽可能小,其中F表示DFT变换。它的核心思想是通过傅里叶变换后,在频域进行乘以Φ,再进行逆傅里叶变换。这又称为频域滤波,因为我们处理的是经过傅里叶变换的信号x而不是x本身。为了找到最佳的对角阵Φ,注意到
令C=Fc和R=Fr为傅里叶变换,同时注意到根据噪声的定义
因此模可以写成
为了确定对角元素
对每个i,解得
但一般情况下我们不知道C或R,所以需要某种近似方法,令X=Fx为傅里叶变换,沿用关于信号与噪声的独立性假设,近似地有
可以把最优参数选为
这里需要用到对噪声等级的先验知识,例如,如果噪声是无关的高斯噪声(独立于信号的正态高斯分布随机数,与信号值相加),可以把
设n为正整数,一维的n阶DCT变换定义为n×n矩阵C,其元素为
其中i,j=0,...,n-1,
定义:设C为如下矩阵,向量x=
注意到C是实数正交矩阵,即转置与逆矩阵相同。
定理(DCT插值定理):令
满足
举例:用DCT变换计算以下点的插值,(0,1),(1,0),(2,-1),(3,0)。
先计算4×4的C矩阵,x=
所以n=4时得到的插值函数为
定理(DCT最小二乘近似定理):设
可以使平方误差
二维DCT变换实际上只是把一维DCT先后应用到二维数据的各个维度上,二维DCT可以用来进行二维网格数据的插值和拟合,计算过程与一维的情形类似。把纵坐标作为第一维,横坐标作为第二维,二维DCT的目的是构造插值函数F(s,t)拟合
二维DCT是将一维DCT变换先后应用到水平和竖直方向上。考虑二维数据
定义:n×n矩阵X的二维DCT变换(2D-DCT)定义为
定义:n×n矩阵Y的二维DCT逆变换是矩阵
正交变换(例如2D-DCT) 与插值有密切的关系,插值的过程是为了恢复原始数据点,先通过DCT变换得到插值系数,然后用这些系数构造插值函数,再在函数上采样得到数据点。由于C是正交阵,
对于
一般对8×8的图像块进行DCT变换,将图像中的能量集中到图像块的左上角。
模q量化
量化:
反量化:
线性量化定义为矩阵
其中0≤k,l≤7,其中常量p称为损失参数。
RGB -> GRAY
RGB -> YUV
不过多赘述
人的听觉系统对频率非常敏感,在压缩和解压过程中产生的任何瑕疵都能被听出来,基于这个原因,音频压缩中通常会用到一些复杂的技巧来掩盖压缩带来的影响。
定义:离散余弦变换版本4(DCT4)是指对n维向量
其中E为n×n矩阵
引理:设
(1)
(2)
我们将用DCT4对应的矩阵E来构造MDCT变换,设n为偶数,我们要用列

定义:令n为正偶数,向量
其中M是n×2n矩阵
其中,0≤i≤n-1,0≤j≤2n-1。
与前面讲到的DCT变换的最大区别是: MDCT把2n维的向量转换为了n维的向量。基于这个原因, MDCT不是可逆的,但是通过重叠这些长度为2n的向量,也可以达到可逆的目的。
我们可以把MDCT的变换矩阵M用DCT4的列来表示,然后用上面的引理来简化,得到
对n=4的MDCT矩阵
令A和B表示DCT4矩阵的左半部分和右半部分,则E=[A|B]。定义一个置换矩阵R,可以把矩阵的各列顺序颠倒:
当一个矩阵右乘以矩阵R,则它的各列的顺序会颠倒。若矩阵左乘以R,则其各行的顺序会颠倒。
利用置换矩阵,可以进一步简化MDCT矩阵的表示
其中AR和BR是列序颠倒后的A和B矩阵。
MDCT变换可以用DCT4来表示,令
为一个2n维的向量,其中
其中E是一个n×n的DCT4变换矩阵,
由于MDCT变换的矩阵M是一个n×2n矩阵,而不是方阵,因此它是不可逆的。但是,两个相邻的MDCT 的总阶数为2n,因此若把它们组合起来,就可以完美地重建出输入数据x。
MDCT的“逆”表示为一个2n×n矩阵
因为E是正交阵,可以得到
在音频压缩中,MDCT变换的向量对应于重叠的数据段。压缩误差会使每段向量之间的连接处产生不连续,而由于向量的长度是固定的,因此这种误差会带来固定频率的噪声。而听觉系统对周期性错误的敏感度要高于视觉系统;实际上,固定周期的错误就是那个频率上的一个声调,而人耳正是感知的各种声调的。如果数据表示采用重叠的方式,令n为正偶数,
为两个2n维的向量,其中每个
把
通过对信号的MDCT进行量化,可以实现音频信号的有损压缩。通过推广图像压缩中使用的量化方法,可以更好地控制信号表示所需的位数。
首先确定一个实数开区间(-L,L),假设我们要用b位来表示(-L,L)上的一个数,而且可以接受一定程度的误差。我们可以用一个位来表示符号,并用b—1位来对数值进行量化。公式为:
(-L,L)的b位量化
量化:
,其中 反量化:
在求解方程中遇到的威尔金森多项式,难以寻找特征值。本节介绍的方法基于对矢量乘上矩阵的高阶幂,这个向量随着幂的升高会变成特征向量,在后面我们将对这个方法进行精化,但是这就是最复杂方法的主要思想。
幂迭代的动机是与矩阵相乘可以将向量推向主特征向量的方向。
定义:令A是一个m×m矩阵,A的占优特征值是量级比矩阵A所有其他特征值都大的特征值λ,如果这样的特征值存在,与λ相关特征向量被称为占优特征向量。
幂迭代本质上是一个在每步上进行归一化的不动点迭代。和FPI一样,该方法线性收敛,意味着在收敛过程中,在每个迭代步骤中误差以常数因子降低。在本节的随后部分,将讨论幂迭代具有二阶收敛的变种,该方法称为瑞利商迭代。
考虑特征值方程xλ=Ax,其中x是特征向量的近似,λ未知,在这种情况下,系数矩阵是n×1矩阵x,法线方程指出最小二乘解是
这就是瑞利商,给定特征向量近似,瑞利商是特征值的最优近似。
幂迭代
给定初始向量
for j = 1,2,3,...
end
定理:令A是一个m×m矩阵,特征值为
幂迭代局限于求解最大(绝对值最大)的特征值,如果幂迭代用于矩阵的逆矩阵,可以找到最小的特征值。
逆向幂迭代
给定初始向量
以及平移s for j = 1,2,3,...
求解
end
为了找出矩阵A在实数s附近的特征值,对
瑞利商可以同逆向幂迭代同时使用,我们知道它会收敛到和s最接近的特征值对应的特征向量,如果这个距离很小则收敛速度很快。如果在这个过程中任何步骤里,近似的特征值已知,可以使用该特征值作为s,以加速收敛。使用瑞利商作为逆向幂迭代中更新的平移会得到瑞利商迭代(RQI)方法。
瑞利商迭代
给定初值
for j=1,2,3,...
求解
end
本节的目标是推导可以一次找出所有特征值的方法,我们从一个用于对称矩阵的方法开始,然后对该方法进行补充,并用于一般矩阵的特征值求解。对称矩阵容易处理是因为它的特征值为实数,其特征向量构成
假设开始有m个两两正交的初始向量,对每个向量使用幂方法一步后得到的向量不再保证彼此正交,所以在每步后重新对m个向量正交化,正交步骤可以看做对结果进行QR分解,如果使用初等基向量作为初始向量,则重新正交化后幂迭代的第一步是
归一化同时迭代(NSI)
设
for j=1,2,3,...
end
在第j步,
xxxxxxxxxxfunction [lam, Q] = nsi(A,k)[m,n] = size(A);Q=eye(m,m);for j = 1:k [Q,R] = qr(A*Q) % QR分解endlam = diag(Q'*A*Q); % 瑞利商
无移动的QR算法
用
xxxxxxxxxxfunction [lam, Q] = unshiftedqr(A,k)[m,n] = size(A);Q=eye(m,m);Qbar = Q; R=A;for j = 1:k [Q,R] = qr(R*Q) % QR分解 Qbar = Qbar * Q;endlam = diag(R*Q); % 瑞利商
定理:假设A是对称m×m矩阵,特征值为
定义:如果矩阵T为上三角矩阵,且在主对角线上可以出现2×2的块,该矩阵具有实数舒尔形式。
定理:令A是实数元素的方阵,则存在正交矩阵Q以及实数舒尔形式的矩阵T,满足
完整的QR算法迭代通过一系列相似变换,将任意矩阵A移动到它的舒尔分解。我们将使用两个步骤,首先将移动的逆向幂迭代思想植入,并加入收缩技术推出移动QR算法,然后推出改进版本允许处理复数特征值。平移版本可以非常直观地写出。每步中都需要使用平移,完成QR分解,然后再平移回来,使用符号表示:
注意到
意味着
原书中有平移QR算法和改进版本的matlab代码
定义:如果
定理:令A是一个方阵,存在正交矩阵Q,满足
可以使用用于构造QR分解的豪斯霍尔德反射构造B,但是在这里有一个主要的差别:现在我们关注在矩阵的左侧和右侧乘上反射子H,这是由于最后想得到一个具有相同的特征值的相似矩阵。正因为如此,我们必须逐步地将0放入矩阵A中。
原书中给出了对应的matlab代码
在
举例:找出矩阵A的奇异值和奇异向量
在对m×n矩阵A分解的过程中,有一个标准的方式来记录所有的信息,生成一个m×m矩阵U,它的列为左奇异向量
引理:令A是一个m×n矩阵,
对于一个m×n的矩阵A,n×n矩阵
如果
如果
注意
)因此是
定理:令A是m×n矩阵,则存在两个正交基:
举例:找出2×2矩阵A的奇异值和奇异向量
的奇异值以降序排列为
求解
注意SVD分解并不唯一,还有另外一种分解
对于对称矩阵,只需计算自身的特征值和特征向量即可。
假设
性质1:矩阵
性质2:如果A为n×n矩阵,则
性质3:如果A为可逆的m×m矩阵,则
性质4:m×n矩阵A可以写成秩为1的矩阵的和
其中r为A的秩,
降维思想是将数据投影到低维空间。假设
SVD则给出了一种直接完成降维的方式,考虑数据向量为m×n矩阵
我们可以将
由于矩阵与
可以利用性质4来压缩矩阵信息,注意到在性质4中的秩1展开中的每一项使用两个向量
如果A是一个实数对称矩阵,SVD可以消减为本章前面讨论的特征值计算问题。在这种情况下单位特征向量构成正交基。如果我们定义矩阵V将单位特征向量保存在列上,则AV=US表示特征向量方程,其中S是对角线矩阵,保存特征值的绝对值。对于U的操作和V一样,但是如果特征值是负,则要变换列的符号,由于U和V是正交矩阵,
是矩阵A的奇异值分解。
对于一个一般的非对称m×n矩阵A,确定SVD有两种不同的方法。最明显的方法是构造
考虑如下矩阵
注意,B是对称的(m+n)×(m+n)矩阵,因而,它具有实数特征值和作为基的特征向量。令[v,w]表示一个(m+n)向量,对应B的特征向量,则
或者Av=λw,左侧乘上
表明w是
因而,另一种更好的计算奇异值和奇异向量的方法是,将对称矩阵B转化为上海森伯格形式。由于对称,上海森伯格形式的矩阵和三对角线矩阵等价,则该方法与平移QR算法一样可以找出特征值,这是奇异值的平方,特征向量中n个最大的项是奇异向量v、尽管这种方式看起来使得矩阵的大小加倍,但是它避免对条件数的不必要提高。
最优化指的是找出实数函数的极大或者极小值,该函数称为目标函数。由于定位函数f(z)的极大值与找出函数-f(x)的极小值等价,在推导计算方式时仅考虑最小化问题就足够了。
定义:当区间[a,b]上只有一个极大或者极小值,并且f在其他点上严格升高或降低,连续函数f(x)被称为区间[a,b]上的单峰函数。
定理:从初始区间[a,b]开始,黄金分割搜索k步之后,最后区间中点和最小值之间的差异为
黄金分割搜索
给出单峰函数f,在区间[a,b]上具有极小值
for i = 1,2,3,...
if f(a+(1-g)(b-a)) < f(a+g(b-a))
b = a+g(b-a)
else
a = a+(1-g)(b-a)
end
end
最终区间[a,b]中包含极小值
持续抛物线插值
从近似极小值r,s,t开始 for i = 1,2,3,...
t = s s = r r = x end
x是对极小值的新的近似。在SPI中,新的x可以替换r,s,t中离当前最远,或者最差的一个点,并重复进行该步骤。对于SPI不能保证收敛,这与黄金分割搜索不同。但是,如果收敛的话,由于该方法更好地使用了函数求值的信息,所以速度会更快。
Nelder-Mead搜索试图将一个多面体滚到一个尽可能低的水平。由于这个原因,它也被称为单纯型下山法。它没有使用目标函数的导数信息。假设需要最小化的函数f具有n个变量,该方法首先需要n+1个属于

具体代码见原书P503页
如果函数连续可微,可以对导数求导,则优化问题可以表示为求解根的问题。
在一个连续可微函数f(x)的局部极小值点
当遇到一个多变量函数时,可以使用多变量的牛顿方法,设置导数为0,有
其中
表示梯度。
使用方程组一章中的牛顿方法可以求解这个问题,设
这是函数f的海森矩阵,因而牛顿步骤为
最速下降,也称为梯度搜索,背后的基本想法是将当前点在在最速下降的方向上移动以找出函数最小值。由于
最速下降
for i = 0,1,2,...
对于标量
最小化
end
当A为对称正定矩阵时,求解Ax=w等价于找出抛物面的极小值,如在二维中,线性方程组
的解是如下抛物面的极小值
这是由于f的梯度可以表示为线性方程组
在极小值处的梯度为0,正定意味着抛物面的凹面向上。
一个关键的观测时线性方程组的余项r=w-Ax为
这意味着
共轭梯度搜索
令
为初始估计,设 for i = 1,2,3,...
使得 最小
end
分块计算A×B,分块的唯一要求是:对A矩阵的各列的分组方法要与B矩阵各行的分组方法一致,如
